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A major challenge of many diffraction calculations, using some form of the Raylcigh- 
Sommcrfcld formulas, is the integration of a highly oscillatory integrand. Here we derive a 
potentially useful alternative form of solution to the Helmholtz equation, which implies a 
sampling theorem for the evaluation of a diffracted scalar field. This alternative solution 
bears close resemblance to the Rayleigh-Sommerfeld diffraction formulas, but instead 
incorporates the boundary conditions digitally. Hence, the integration is replaced by a 
simple summation. This formulation may be more efficient for accurate computer-based 
calculation of the diffracted scalar field. © 2013 Optical Society of America 



1. Introduction 

The diffraction of monochromatic electromagnetic radiation by obstacles much larger than the wavelength is 
often treated satisfactorily by the Scalar Diffraction Theory. That is, one solves the Helmholtz equation with 
certain boundary conditions. We let U be any component of either the electric or magnetic fields, defined 
in the half-space z > 0, free of charges or currents. If the scalar field U is simple harmonic in time with 
frequency u, then it satisfies the Helmholtz equation, 

V 2 U + k 2 U = (1) 

where V 2 is the Laplacian differential operator, k = ut/c = 2ir/\ is the wavenumber, c is the speed of light 
and A is the wavelength. The boundary conditions are usually given as follows [T] . For some known function 
Uo(x, y), the scalar field U on the boundary plane z — is U(x, y, 0) = Uo(x, y). If r is the distance from the 
origin to any point in the half-space, we assume the Sommcrfcld radiation condition, lim r _ ! . 0O U = 0(l/r), 
such that the total energy is bounded. Additionally, we assume that all radiation is propagating in the 
positive z-direction. 

One solution to this boundary value problem is the Rayleigh-Sommerfeld diffraction formula of the first 
kind, which can be written as [T], 

U(x,y,z) = =± jju {x',y')^^-dx'dy' (2) 

where R is the distance between the points P and Q with coordinates (x, y, z) and (x', y', 0), respectively, and 
the integration is taken over the entire plane z = 0. This integral can be evaluated exactly only for a small 
class of boundary value functions Uo(x,y). Under certain additional constraints, it can be approximated 
by the Fresnel or Fraunhofer diffraction integrals [2H4], which essentially reduces the problem to Fourier 
analysis. However, the diffracted field is often found from Eq. ([2]), given some general incident field Uo(x,y), 
by evaluating the integral numerically using a computer. 
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The principal challenge of evaluating the right side of Eq. numerically is creating a sampling mesh 
fine enough to accurately represent the highly oscillatory integrand. The exponential term oscillates spatially 
with a mimimum wavelength of A. Thus, it seems necessary that the spacing between mesh points be much 
smaller than A. Here we show that a sampling mesh with spacing A/2 between points is sufficient to obtain 
high accuracy for most practical applications. Moreover, if the field is composed only of homogeneous plane 
waves, then this produces the exact result. 

We propose a modified version of the Rayleigh-Sommerfeld solution, which incorporates the boundary 
conditions discretely. To arrive at this formulation, we employ the angular spectrum method [5] of field 
analysis. Both the two-dimensional and three-dimensional versions of the result are discussed. 

2. The Angular Spectrum: Homogeneous and Inhomogeneous Wave Components 

The Rayleigh-Sommerfeld solution to the Helmholtz equation is obtained by using the Kirchhoff integral 
formula in conjunction with the Green's function solution pQ. Yet, there is an alternative method for solving 
Eq. (1), referred to as the angular spectrum method. It can be shown [TJ[S] that the general solution to Eq. 
(2) in the half-space z > with the given boundary conditions can be written as 

U(x,y,z) = // A(u,v) exp(ik[ux + vy + wz]) du dv (3) 



where w = +yl — u 2 ~ v 2 , and 

A(u,v) = \~ 2 I I Uq(x, y) exp(— ik[ux + vy\) dx dy (4) 



is the so-called angular spectrum of the field U, as a function of the direction cosines u and v. Evidently, the 
angular spectrum function is just the Fourier transform of the boundary value function Uo(x,y). For each 
plane wave component, 

Uu,v = A(u, v) exp(ik[ux + vy + wz]) 

in the integrand of Eq. ((5J), there are two physically distinct possibilities: u 2 + v 2 < 1 and u 2 + v 2 > 1. In 
the former, w is real and the plane wave is called homogeneous. In the latter case, w is imaginary and the 
plane wave is called inhomogeneous, 

exp(ik[ux + vy + wz]) = cxp(—k\w\z) exp(ik[ux + vy]). (5) 

This latter kind of wave is also known as an evanescent wave, because it propagates in the plane z = and 
vanishes rapidly in the positive z direction. Based on this consideration, wc can express the field U as the 
sum of its homogeneous and inhomogeneous plane wave components, 

U(x,y,z) = U H (x,y,z) + U I (x,y,z) (6) 

where the homogeneous component is 

Uh(x,u,z)= // A(u,v) exp(ik[ux + vy + wz]) du dv, (7) 

J Ju 2 + v 2 <l 

and the inhomogeneous component is 

Ui(x,y,z)= // A(u,v) exp(—k\w\z) exp(ik[ux + vy]) du dv. (8) 

J Ju 2 +v 2 >l 
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Let us place an upper bound on \Ui\ for any (real) value of kz. Applying Schwarz' inequality to the right 
side of Eq. (|5J) gives 

A(u, v) exp(—k\w\z) exp(ik[ux + vy]) du dv\ 

U 2 +V 2 >1 ^ 

< \A(u,v)\ 2 du dv ■ // exp(-2k\w\z) du dv 

JJu 2 +v 2 >l JJu 2 +v 2 >l 

The second integral on the right side is easily evaluated. Switching from Cartesian to polar coordinates in 
the space of the direction cosines, with radial coordinate p = \Ju? + v 2 , this integral takes the form, 

(10) 



(11) 



/ exp(— 2fc|w|z) du dv = 2ir / exp(— 2kz\J p 2 — 1) p dp 

'u 2 +v 2 >l Jl 

Then, using the substitution t 2 = p 2 — 1 yields 

/>oo poo 

2tt / exp(— 2kz\/ p 2 — 1) p dp = 2tt / cxp(— 2kzt) t dt 
Ji Jo 

Finally, applying integration by parts to the right side has the result, 



lu 2 +v 2 >l 

Combining this with Eq. (10), we obtain 



exp(—2k\w\z) du dv = (12) 
2k z z z 



\Ui(x, y ,z)\<± ?L [[ \A(u,v)\ 2 du dv. (13) 

kz V 2 JJ U 2 +V 2 >1 

In other words, for any given angular spectrum function A(u,v), the optical field is given asymptically by 

U(x,y,z) = U H (x,y,z) + 0(l/kz). (14) 

as kz — > oo. In most applications, the value of kz is very large, e.g. ~ 10 6 for visible light. It follows that the 
inhomogeneous plane wave component E/j is usually nearly zero in almost all of the half-space z > 0. This 
contribution will be neglected in the following analysis. 

In dealing only with the homogenous plane wave component, it is desirable to express Ujj in terms of 
its values at the boundary z = 0. The function Uh (x, y, 0) can be expressed in simple terms of the known 
function Uq(x, y). It follows from elementary Fourier analysis and the convolution theorem that 



U H (x,y,0) = U (x,y) * Jl(fc /f + f } (15) 

Ay x £ + y z 

where * is the convolution operator and J\ is the first order Bessel function of the first kind. It is easy to see 
that as the wavelength A approaches zero, the second term in the convolution becomes the two dimensional 
Dirac delta function 6 2 (x,y), with the property Uq * S 2 = Uq. In regions where Uo(x,y) varies with length 
scales much larger than A, this implies Uh(x, y, 0) ~ Uo(x, y). This approximation will fail at the boundaries 
of the obstructing obstacle where the field has a discontinuity. However, if the size of the obstacle is much 
larger than the wavelength, then the field near this boundary will contribute negligibly to the diffracted field 
at points far away from the obstacle. 
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3. Diffraction in Two Dimensions 

Many configurations of diffraction are two-dimensional in form, e.g. diffraction by an infinitely long slit. Then 
as a first case, let us assume that the boundary field Uq(x, y) = Uq(x) on the plane z = has no dependence 
on the y-coordinatc. Then, from Eq. (4) the angular spectrum can be written as 

A(u,v) = S(v)A'(u) (16) 

where 5 is the (one dimensional) Dirac delta function. Inserting this into Eq. (7) then yields 

Uh(x,V, z) = Uh(x, z) = I A(u) exp(ik[ux + wz]) du (17) 



where here w = +y/l — u 2 and the dash is dropped from the angular spectrum function A(u). As expected, 
the homogeneous component in z > also has no y-dependence. The same argument obviously applies to 
Ul 

To evaluate the integral in Eq. (flTj) . the function A(u) need only be defined on the interval u € [—1,1]. 
Within this domain, let us expand A(u) as the Fourier series, 

oo 

A(u) = c m exp(— iirmu) (18) 

m— — oo 

which is always possible for any sufficiently well-behaved function A(u) [6]. The expansion coefficients c m 
are found by the inversion formula, 



1 ^ 



A(u) exp(i7T77ra) du. (19) 



-l 



Comparing Eq. (|19j) to Eq. ()17|) . it is immediately apparent that 



W^,0). (20) 



„. 2 V 2 

Combining this relation with Eqs. (16) and (17) therefore has the result 

■sr-> ( TnX \ „ / mX \ 
U H (x,z)= J2 Uh ( — > ) G ( x -—> z ) ( 21 ) 

m— — oo \ / \ / 

where 

1 f 1 

G{x,z) = — I cxp(ik[ux + wz]) du. (22) 

2 J-i 

We can compare Eq. (|2~Tj) to Eq. (2). For the homogeneous component Uji of the optical field U, we replace 
the integral of Eq. ([2]) with the simple sum in Eq. (|21|). That is, instead of requiring knowledge of the whole 
function Uh on the boundary z — 0, we only need to know its values at uniformly spaced points x = mX/2 
for all integer m. Moreover, if the field is only non-zero for a finite region on z = 0, then one only needs a 
finite set of known values of Uh on z = to perfectly compute Uh anywhere else in the half-space z > 0. 

Let us now evaluate the integral on the right side of Eq. (f22j) for a clearer expression of the kernel function 
G(x, z). Firstly, consider the two important limiting cases: z = and z » X. The first case is simple, 

1 f 1 

G(x,0) = — J cxp(ikux) du = sin(kx)/kx, (23) 
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which reduces the right side of Eq. (j!?Tj) to Whittaker's cardinal series expansion [7HTU] of the function 
Uh(x,0). For the second case (z >> A, or equivalently kz » 1), we make the trigonometric substitution, 
u = sin 77, such that w = cos 77 and using basic trigonometric identities Eq. (|22l) becomes 

1 r /2 

G(x,z) = — / cxp(ifcr cos[?7 — 8]) cos 77 dr/, (24) 

2 J-tt/2 

switching to polar coordinates: x = rsu\8 and y = rcosd. The restriction z > implies that the angle 8 is 
confined to 6* e [— 7r/2, 7r/2]. Because we have kz >> 1, the exponential term oscillates rapidly, except near 
the stationary points of cos[rj — 8]. These stationary points are 77 = 8 + lir for any integer I. However, within 
the integration interval [— tt/2, 7t/2] there can only be one stationary point, namely 77 = 9. Furthermore, near 
this point cos[?7 — 8] ps 1 — [77 — 8] 2 /2. Invoking the principle of stationary phase [TT], it follows that 



n/2 r , exp(7;fcr) 

exp(ifcr cos[77 — wj) cos 77 (177 sa cos( 

■tt/2 2 



exp 



^77. 



This well known integral is easily evaluated, resulting in the approximation 

<tt exp(ifcr) 



G(x, z » A) 



'ikr 



■ cos 8. 



(25) 



(26) 



This function is just the asymptotic form (for large kr) of the two dimensional field established by an 
oscillating dipole at the origin, with the dipole moment aligned with the x-axis. Combining Eqs. (|21[) and 
(l26l). we obtain 



U H (x,z » A) 




mX \ cxp(tkR m ) 
~^~'°) 7^= — cos ^ 

* / V -Km 



(27) 



where i? r „ is the distance between the points at (x,z) and (mA/2,0), and 8 m is the angle between the ray 
connecting these two points and the z-axis as illustrated in Fig. 1. 



(x, z) 




z = 



Fig. 1. A two-dimensional electromagnetic disturbance, composed solely of homogeneous 
plane waves and satisfying the given boundary conditions, is uniquely determined at any 
point (x, z) in the half-space z > by its values on the boundary z = at points separated 
by A/2. The field is constructed by placing individual oscillating dipoles at each of these 
boundary points and adding up all contributions at the point of interest. 



To the knowledge of the author, there is no exact closed-form expression for the kernel function G(x, z) 
in terms of elementary functions for the two-dimensional case. 
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4. Diffraction in Three Dimensions 



Wc now examine the more general case of diffraction in three dimensions, dealing only with the homogeneous 
component Uh of the field U. Starting from Eq. (7), the evaluation of the integral only requires that the 
angular spectrum function A(u, v) be defined on the unit disc u 2 + v 2 < 1. Let the function f(u, v) be defined 
on the square domain, in which both u and v are bound to the interval [—1, 1], such that 



/<„,„) =< »:+«:;; > <*» 

u + v > 1 



Within this square domain, the function f(u,v) can be expanded as the two-dimensional Fourier series, 



oo 



f(u,v)= } c m ,„ exp(-i7r[mw + nv]). (29) 



m,n— — oo 



As before, the expansion coefficients are 

c m .n = \ J_i /_i /("", w) exp(i7r[mu + nv]) du dv 

= \ IJ u3+«2< 1 ^( W : u ) exp(wr[?mt + ra]) d« (30) 

- i U H { — , — > U J 

where the third equality follows from the second upon comparing with Eq. (|7|)- Because f(u,v) = A(u,v) 
within the unit disc, the function A(u,v) can be replaced by f(u,v) in Eq. ([7]). Combining this with Eqs. 
([29]) and (J3DJ) then yields 

U H ( X ,y,z)= f uJ^^GL-^y-^z) (31) 



m,n— — oo 



where 

G(x,y,z) = — // exp(ifc[wa; + + mz]) c?w dw. (32) 

4 J Ju 2 +v 2 <l 

is the three dimensional kernel function. 

To evaluate the integral in Eq. (|3"2"j) . we make use of the coordinate transformations illustrated by Fig. 2 
and described as follows. The approach is similar to that in Ref. |12j . but here a more careful transformation 
of the limits of integration is required. We first make the transformation [x, y, z) — > (x', y' , z'), consisting of 
a rotation about the origin, such that the point P with coordinates (x, y, z) lies on the z'-axis in the new 
coordinate system. That is, x' = y' = for point P. The same transformation is also applied to the direction 
cosines, (u, v, w) — > (u',v', w'). This change of variables preserves the dot product, 

ux + vy + wz = u x + v y + w z = w'r (33) 



where r = y x 2 + y 2 + z 2 . The ratio of the differential areas du dv and du' dv' (given by the determinant of 
the Jacobian) is cos#, where 8 = cos _1 (z/r) is the angle between point P and the z-axis. With this, 

\ J J e M ik[ux + vy + wz] ) dudv^ J J exp (i Wr) dv' . (34) 

This form suggests a switch to spherical coordinates in the space of the direction cosines, 

u' = sin n cosip 

v' = sinr; sin^ (35) 
w' = cos r\ 



ti 



where rj and ip are the azimuthal and polar angles, respectively, of the point with coordinates (it', v', w'). 
The differential area transforms once more as 



du' dv' = sin 77 cos 77 drj dip. 



(36) 



The limits of integration, defined by ti 2 + jj 2 < 1 in the initial coordinate system, transform accordingly. 
The polar angle ip is integrated from to 2tt. The azimuthal angle r\ is integrated from to ir/2 + A77, 
where Arj is the angle between the planes z = and z' = as seen from the origin at angle ip. Applying 
the spherical law of sines [12] to the geometry shown in the inset of Fig. 2, it can be easily verified that 
A77 = sin (cos ^ sin 0). Therefore, the integral in Eq. (|52")l becomes 



3 / fu 2 +v 2 <i exp(ifc[wx + vy + wz}) du dv 



TfH^ jt^^^eMikrcosrismrj dr, d^ 



1 d r27T p-^+sin 1 (cos t/j sin 8) 
4ik7z JO JO 



exp(ifcrcos?7) sin 77 drj dip. 



(37) 

The chain rule is applied, cos# (d/dr) = (dr /dz)(d/dr) = d/dz, to obtain the second equality. The integral 




Fig. 2. Geometry of the coordinate transformation. The point P in the new coordinate 
system lies on the z'-axis. The point given by coordinates (u,v,w) lies on the unit half- 
sphere w > 0. The shaded region is the plane z = 0. In the inset, right angles of the 
spherical triangle are denoted by small squares. 
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(38) 



on the right side of Eq. (37) can now be evaluated with the substitution r = cos 77, yielding 

1 ttv,*+v*<i «xp(ik[wc + vy + wz\) du dv = J- sin ecos^ exp(iArr) dr dtp 

1 d r27r cxp(ifcr)— cxp( — ikr sin 9 cos i i 
~ Welz JO Tkr a ^ 

and finally, 

1 /A\ 2 d exp(ifcr) - J (fcrsin6») 
G( W ) = --(^-j . (39) 

where Jo is the zeroth order Bessel function of the first find, which admits the identity [14] , 

1 f 27T 

Jo(s) = — / exp(zscos0) d(j>, (40) 
27r Jo 



used to obtain Eq. (|39|) from Eq. (|38|) . Note that, whereas we were only able to obtain an approximate 
simplification of the kernel function G in the two-dimensional case, Eq. (f3U)) is an exact expression of this 
function for three dimensions. Applying this to Eq. (|3 1[) gives the discretized diffraction formula, 



U H ( X ,y,z) = --[-) ^ U H [—,—,0)- — (41) 



m,n— — oo 



where R m ,n — y (x — mA/2) 2 + (y — nX/2) 2 + z 2 is the distance between the point at (mA/2, nA/2, 0) and 
point P, and 6 m . n = cos~ 1 (z/i? I „ ! „) is the angle between the ray connecting these two points and the z-axis, 
as illustrated in Fig. 3. Equation (|4ip bears much resemblance to the Raylcigh-Sommcrfcld integral formula 
of Eq. ([2]), with the exception of the Bessel function term, and more importantly the integral is replaced by 
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(x,y,z) 




Fig. 3. The homogeneous component Uh of the diffracted field at the point (x, y, z) in the 
half-space z > 1 is completely determined by its values on the rectangular lattice of points 
at (mA/2, nA/2, 0), where m and n are integers. It is found by placing secondary sources at 
each lattice point and adding up each contribution at the point of interest. 

5. Discussion and Conclusions 

We present a solution to the scalar Hclmholtz equation, the cornerstone of Scalar Diffraction Theory, which 
inherently implies a potentially useful sampling theorem regarding the propagation of electromagnetic ra- 
diation in free space. Summarily, if the scalar field is composed only of homogeneous plane waves, then 
it can be perfectly reconstructed from the knowledge of its values on a lattice of points on the boundary 
z = 0. The lattice of points has uniform spacing of A/2 between collinear points. If the field is composed 
of inhomogeneous plane waves, then one can always choose some point z such that this contribution to the 
total field becomes negligibly small. 

Going further, the solution presented in Eq. pTj) can be straightforwardly interpreted as a matrix equa- 
tion, if the field were to be evaluated on a secondary lattice of points for some z > 0. Posing the diffraction 
calculation in this form would be well-suited for numerical algorithms performed by a computer. The opera- 
tion could possibly be optimized for efficiency, analogous to the Fast Fourier Transform. Such a development 
would significantly reduce the time needed to calculate the diffracted field in many applications. 
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